library(foreign)
library(ggplot2)
library(colorspace)

# setwd("")


summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
                      conf.interval=.95, .drop=TRUE) {
  require(plyr)
  
  # New version of length which can handle NA's: if na.rm==T, don't count them
  length2 <- function (x, na.rm=FALSE) {
    if (na.rm) sum(!is.na(x))
    else       length(x)
  }
  
  # This does the summary. For each group's data frame, return a vector with
  # N, mean, and sd
  datac <- plyr::ddply(data, groupvars, .drop=.drop,
                       .fun = function(xx, col) {
                         c(N    = length2(xx[[col]], na.rm=na.rm),
                           mean = mean   (xx[[col]], na.rm=na.rm),
                           sd   = sd     (xx[[col]], na.rm=na.rm)
                         )
                       },
                       measurevar
  )
  
  # Rename the "mean" column    
  datac <- plyr::rename(datac, c("mean" = measurevar))
  
  datac$se <- datac$sd / sqrt(datac$N)  # Calculate standard error of the mean
  
  # Confidence interval multiplier for standard error
  # Calculate t-statistic for confidence interval: 
  # e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
  ciMult <- qt(conf.interval/2 + .5, datac$N-1)
  datac$ci <- datac$se * ciMult
  
  return(datac)
}


swiss <- read.csv("SwissReplicationData.csv",encoding="UTF-8")


#### FIGURE 3: KNOWLEDGE ABOUT SWISS SCHEME AMONG DECEMBER 2019 SURVEY RESPONDENTS #### 

ggplot(swiss, aes(x=taxpresent, fill=taxpresent)) + 
  geom_bar(aes(y = (..count..)/sum(..count..)),show.legend=FALSE)+
  scale_y_continuous(limits=c(0,1)) +
  scale_x_discrete(labels=c("yes","no","don't know"))+
  scale_fill_manual(values=c("#009E73", "#000000","#000000"))+
  ylab("Fraction of survey respondents") +
  xlab("Belief that Switzerland has a carbon tax on fossil fuels") +
  theme_bw()

dev.copy2pdf(file="PaperFiguresFinal/Figure3_topleft.pdf")


ggplot(swiss, aes(x=taxuse, fill=taxuse)) + 
  geom_bar(aes(y = (..count..)/sum(..count..)),show.legend=FALSE)+
  scale_y_continuous(limits=c(0,1)) +
  scale_x_discrete(labels=c("public rebates","renewable energy",
                                              "gov. spending.","don't know"))+
  scale_fill_manual(values=c("#009E73", "#000000","#000000", "#000000"))+
  ylab("Fraction of survey respondents") +
  xlab("Most carbon tax revenue is directed towards") +
  theme_bw()

dev.copy2pdf(file="PaperFiguresFinal/Figure3_topright")


ggplot(swiss[-which(is.na(swiss$rebatepurpose)),], aes(x=rebatepurpose, fill=rebatepurpose)) + 
  geom_bar(aes(y = (..count..)/sum(..count..)),show.legend=FALSE)+
  scale_y_continuous(limits=c(0,1)) +
  scale_x_discrete(labels=c("cantonal tax", "national tax", "health bill", "pension contr.", "electr. bills", "don't know"))+
  scale_fill_manual(values=c("#000000","#000000","#009E73", "#000000","#000000", "#000000"))+
  ylab("Fraction of survey respondents") +
  xlab("Tax revenue redistributed to reduce...") +
  theme_bw()

  dev.copy2pdf(file="PaperFiguresFinal/Figure3_bottomleft")
  

ggplot(swiss[-which(is.na(swiss$rebateguess)),], aes(x=rebateguess, fill=rebateguess)) + 
  geom_bar(aes(y = (..count..)/sum(..count..)),show.legend=FALSE)+
  scale_y_continuous(limits=c(0,1)) +
  scale_x_discrete(labels=c("None", "less than 3CHF", "3-10CHF", "10-15CHF", "Don't know"))+
  scale_fill_manual(values=c("#000000","#000000","#009E73", "#000000", "#000000"))+
  ylab("Fraction of survey respondents") +
  xlab("Perceived monthly rebate size") +
  theme_bw()

dev.copy2pdf(file="PaperFiguresFinal/Figure3_bottomright")



#### FIGURE 5: RESULTS OF SWISS SURVEY EXPERIMENT #### 

swiss_cpsupport <- summarySE(data=swiss, measurevar=c("support"), groupvars=c("exptreat"), na.rm=TRUE)
swiss_cpsupport <- swiss_cpsupport[-3,]

swiss_cpsupportleft <- summarySE(data=swiss[which(swiss$partyideo=="left party"),], measurevar=c("support"), groupvars=c("exptreat"), na.rm=TRUE)
swiss_cpsupportleft <- swiss_cpsupportleft[-3,]

swiss_cpsupportcentre <- summarySE(data=swiss[which(swiss$partyideo=="centre party"),], measurevar=c("support"), groupvars=c("exptreat"), na.rm=TRUE)
swiss_cpsupportcentre <- swiss_cpsupportcentre[-3,]

swiss_cpsupportright <- summarySE(data=swiss[which(swiss$partyideo=="right party"),], measurevar=c("support"), groupvars=c("exptreat"), na.rm=TRUE)
swiss_cpsupportright <- swiss_cpsupportright[-3,]

swiss_cpsupport_graph <- rbind (swiss_cpsupport, swiss_cpsupportleft, swiss_cpsupportcentre, swiss_cpsupportright)
swiss_cpsupport_graph$group <- c("all", "all", "left", "left", "centre", "centre", "right", "right")
swiss_cpsupport_graph$group <- ordered(swiss_cpsupport_graph$group, levels=c("all","right", "centre", "left"))

quartz(width=7, height=5.5)

ggplot(swiss_cpsupport_graph, aes(x=as.factor(exptreat), y=support, color=group)) + 
  geom_point(position=position_dodge(.4), size=4)+
  scale_y_continuous(limits=c(1,4)) +
  scale_color_manual(values=c("#000000","#56B4E9","#009E73", "#D55E00"))+
  ylab("Support for existing policy (4 point scale)") +
  xlab("Treatment Status") +
  theme_bw()+
  geom_errorbar(aes(ymin=support-ci, ymax=support+ci),
                width=.2,
                size=1,
                position=position_dodge(.4))

dev.copy2pdf(file="PaperFiguresFinal/Figure5_top.pdf")



## low only

swiss_cpincreaseLOW <- summarySE(data=swiss[which(swiss$hl.dummy==0),], measurevar=c("increase"), groupvars=c("exptreat"), na.rm=TRUE)
swiss_cpincreaseLOW <- swiss_cpincreaseLOW[-3,]

swiss_cpincreaseLOWleft <- summarySE(data=swiss[which(swiss$partyideo=="left party"&swiss$hl.dummy==0),], measurevar=c("increase"), groupvars=c("exptreat"), na.rm=TRUE)
swiss_cpincreaseLOWleft <- swiss_cpincreaseLOWleft[-3,]

swiss_cpincreaseLOWcentre <- summarySE(data=swiss[which(swiss$partyideo=="centre party"&swiss$hl.dummy==0),], measurevar=c("increase"), groupvars=c("exptreat"), na.rm=TRUE)
swiss_cpincreaseLOWcentre <- swiss_cpincreaseLOWcentre[-3,]

swiss_cpincreaseLOWright <- summarySE(data=swiss[which(swiss$partyideo=="right party"&swiss$hl.dummy==0),], measurevar=c("increase"), groupvars=c("exptreat"), na.rm=TRUE)
swiss_cpincreaseLOWright <- swiss_cpincreaseLOWright[-3,]

swiss_cpincrease_graphLOW <- rbind (swiss_cpincreaseLOW, swiss_cpincreaseLOWleft, swiss_cpincreaseLOWcentre, swiss_cpincreaseLOWright)
swiss_cpincrease_graphLOW$group <- c("all", "all", "left", "left", "centre", "centre", "right", "right")
swiss_cpincrease_graphLOW$group <- ordered(swiss_cpincrease_graphLOW$group, levels=c("all","right", "centre", "left"))



quartz(width=7, height=5.5)

ggplot(swiss_cpincrease_graphLOW, aes(x=as.factor(exptreat), y=increase, color=group)) + 
  geom_point(position=position_dodge(.4), size=4)+
  scale_y_continuous(limits=c(1,4)) +
  scale_color_manual(values=c("#000000","#56B4E9","#009E73", "#D55E00"))+
  ylab("Support for small carbon tax increase (4 point scale)") +
  xlab("Treatment Status") +
  theme_bw()+
  geom_errorbar(aes(ymin=increase-ci, ymax=increase+ci),
                width=.2, 
                size=1,
                position=position_dodge(.4))


dev.copy2pdf(file="PaperFiguresFinal/Figure5_bottomleft.pdf")



## high only: 

swiss_cpincreaseHIGH <- summarySE(data=swiss[which(swiss$hl.dummy==1),], measurevar=c("increase"), groupvars=c("exptreat"), na.rm=TRUE)
swiss_cpincreaseHIGH <- swiss_cpincreaseHIGH[-3,]

swiss_cpincreaseHIGHleft <- summarySE(data=swiss[which(swiss$partyideo=="left party"&swiss$hl.dummy==1),], measurevar=c("increase"), groupvars=c("exptreat"), na.rm=TRUE)
swiss_cpincreaseHIGHleft <- swiss_cpincreaseHIGHleft[-3,]

swiss_cpincreaseHIGHcentre <- summarySE(data=swiss[which(swiss$partyideo=="centre party"&swiss$hl.dummy==1),], measurevar=c("increase"), groupvars=c("exptreat"), na.rm=TRUE)
swiss_cpincreaseHIGHcentre <- swiss_cpincreaseHIGHcentre[-3,]

swiss_cpincreaseHIGHright <- summarySE(data=swiss[which(swiss$partyideo=="right party"&swiss$hl.dummy==1),], measurevar=c("increase"), groupvars=c("exptreat"), na.rm=TRUE)
swiss_cpincreaseHIGHright <- swiss_cpincreaseHIGHright[-3,]

swiss_cpincrease_graphHIGH <- rbind (swiss_cpincreaseHIGH, swiss_cpincreaseHIGHleft, swiss_cpincreaseHIGHcentre, swiss_cpincreaseHIGHright)
swiss_cpincrease_graphHIGH$group <- c("all", "all", "left", "left", "centre", "centre", "right", "right")
swiss_cpincrease_graphHIGH$group <- ordered(swiss_cpincrease_graphHIGH$group, levels=c("all","right", "centre", "left"))


quartz(width=7, height=5.5)

ggplot(swiss_cpincrease_graphHIGH, aes(x=as.factor(exptreat), y=increase, color=group)) + 
  geom_point(position=position_dodge(.4), size=4)+
  scale_y_continuous(limits=c(1,4)) +
  scale_color_manual(values=c("#000000","#56B4E9","#009E73", "#D55E00"))+
  ylab("Support for large carbon tax increase (4 point scale)") +
  xlab("Treatment Status") +
  theme_bw()+
  geom_errorbar(aes(ymin=increase-ci, ymax=increase+ci),
                width=.2, 
                size=1,
                position=position_dodge(.4))


dev.copy2pdf(file="PaperFiguresFinal/Figure5_bottomright")


# significance tests
mean(swiss$support[which(swiss$partyideo=="left party")], na.rm=TRUE)
mean(swiss$support[which(swiss$partyideo=="right party")], na.rm=TRUE)

summary(lm(support~exptreat, data=swiss))
summary(lm(support~exptreat, data=swiss[which(swiss$partyideo=="left party"),]))
summary(lm(support~exptreat, data=swiss[which(swiss$partyideo=="centre party"),]))
summary(lm(support~exptreat, data=swiss[which(swiss$partyideo=="right party"),]))

summary(lm(increase~exptreat, data=swiss))
summary(lm(increase~exptreat, data=swiss[which(swiss$partyideo=="left party"),]))
summary(lm(increase~exptreat, data=swiss[which(swiss$partyideo=="centre party"),]))
summary(lm(increase~exptreat, data=swiss[which(swiss$partyideo=="right party"),]))

summary(lm(increase~exptreat, data=swiss[which(swiss$hl.dummy==0),]))
summary(lm(increase~exptreat, data=swiss[which(swiss$hl.dummy==0&swiss$partyideo=="left party"),]))
summary(lm(increase~exptreat, data=swiss[which(swiss$hl.dummy==0&swiss$partyideo=="centre party"),]))
summary(lm(increase~exptreat, data=swiss[which(swiss$hl.dummy==0&swiss$partyideo=="right party"),]))

summary(lm(increase~exptreat, data=swiss[which(swiss$hl.dummy==1),]))
summary(lm(increase~exptreat, data=swiss[which(swiss$hl.dummy==1&swiss$partyideo=="left party"),]))
summary(lm(increase~exptreat, data=swiss[which(swiss$hl.dummy==1&swiss$partyideo=="centre party"),]))
summary(lm(increase~exptreat, data=swiss[which(swiss$hl.dummy==1&swiss$partyideo=="right party"),]))





#### SI SECTION 16: EXPERIMENTAL BALANCE #### 

library(coin)
library(arsenal)
my_labels <- list(
  treat1 = "Treatment",
  gender = "Gender", 
  age = "Age",
  education = "Education",
  knowledgeCT = "Prior knowledge about CO2 levy",
  climscep = "Belief in antropogenic climate change"
)

table_treat <- tableby(exptreat ~ 
                         gender + age + education +
                         knowledgeCT + climscep, 
                       data = swiss
)

summary(table_treat, 
        labelTranslations = my_labels, text = "latex",
        title = "Descriptive Statistics")



#DIM-Tests

#generate Dummies per category
swiss$male <- NA
swiss$male[swiss$gender=="Mann"]<- 1
swiss$male[swiss$gender=="Frau"]<- 0

swiss$female <- NA
swiss$female[swiss$gender=="Mann"]<- 0
swiss$female[swiss$gender=="Frau"]<- 1


swiss$age1 <- ifelse(swiss$age == '18 bis 24 Jahre', 1, 0)
swiss$age2 <- ifelse(swiss$age == '25 bis 34 Jahre', 1, 0)
swiss$age3 <- ifelse(swiss$age == '35 bis 44 Jahre', 1, 0)
swiss$age4 <- ifelse(swiss$age == '45 bis 54 Jahre', 1, 0)
swiss$age5 <- ifelse(swiss$age == '55 bis 64 Jahre', 1, 0)
swiss$age6 <- ifelse(swiss$age == '65 bis 74 Jahre', 1, 0)
swiss$age7 <- ifelse(swiss$age == '75+ Jahre', 1, 0)

swiss$edu1 <- ifelse(swiss$education == 'Sec. I', 1, 0)
swiss$edu2 <- ifelse(swiss$education == 'Sec. II', 1, 0)
swiss$edu3 <- ifelse(swiss$education == 'Tertiary', 1, 0)

swiss$know1 <- ifelse(swiss$knowledgeCT == 'Does not know', 1, 0)
swiss$know2 <- ifelse(swiss$knowledgeCT == 'Knows little', 1, 0)
swiss$know3 <- ifelse(swiss$knowledgeCT == 'Knows', 1, 0)


#Test mean difference between treatment and control group

treat <- subset(swiss, exptreat=="treat")
control <- subset(swiss, exptreat=="control")

t.test(treat$male,control$male)

t.test(treat$age1,control$age1)
t.test(treat$age2,control$age2)
t.test(treat$age3,control$age3)
t.test(treat$age4,control$age4)
t.test(treat$age5,control$age5)
t.test(treat$age6,control$age6)
t.test(treat$age7,control$age7)

t.test(treat$edu1,control$edu1)
t.test(treat$edu2,control$edu2)
t.test(treat$edu3,control$edu3)

t.test(treat$know1,control$know1)
t.test(treat$know2,control$know2)
t.test(treat$know3,control$know3)

t.test(treat$climscep,control$climscep)



#### SI SECTION 25: SAMPLE REPRESENTATIVENESS #### 

prop.table(table(swiss$gender))
prop.table(table(swiss$age))
prop.table(table(swiss$canton))
prop.table(table(swiss$education))


#### SI SECTION 27: VARIABLES AND DESCRIPTIVE STATISTICS #### 

my_controls <- tableby.control(
  test = T,
  total = T,
  numeric.test = "kwt", cat.test = "chisq",
  numeric.stats = c("meansd", "medianq1q3", "range", "Nmiss2"),
  cat.stats = c("countpct", "Nmiss2"),
  stats.labels = list(
    meansd = "Mean (SD)",
    medianq1q3 = "Median (Q1, Q3)",
    range = "Min - Max",
    Nmiss2 = "Missing"
  )
)

my_labels <- list(
  exptreat = "Treatment 1",
  know_ct_engl = "Knowledge about the existence a CO2 levy",
  use_ct_engl = "Knowledge about revenue recycling",
  knowledge3 = "Knowledge about how rebate is redistributed",
  knowledge4 = "Knowledge about size of rebate",
  knowledgeCT = "Prior knowledge (Index)",
  gender = "Gender", 
  age = "Age",
  canton = "Home canton",
  education = "Education",
  inc_cat = "Income",
  climscep = "Belief in antropogenic climate change"
)


table_descr <- tableby( ~  exptreat + 
                          know_ct_engl + use_ct_engl + knowledge3 + knowledge4 +
                          knowledgeCT + gender + age + canton + education +
                          inc_cat +  climscep, 
                        data = swiss,
                        control = my_controls)


summary(table_descr, 
        labelTranslations = my_labels, text = "latex",
        title = "Descriptive Statistics")

































